#ifndef RSIM_EIG_H
#define RSIM_EIG_H

/*-------------------------------------------------------------------------- *
* File: Eig.h                                                                *
* Authors: Kishor Bhalerao                                                   *
* Email : kishor8dm@gmail.com                                                *
* Contributors:                                                              *
* Email:                                                                     *
* copyright (c) 2010 Authors.                                                *
*                                                                            *
* This program is free software: you can redistribute it and/or modify it    *
* under the terms of the GNU General Public License as published by the Free *
* Software Foundation, either version 3 of the License, or any later version.*
*                                                                            *
* This program is distributed in the hope that it will be useful, but        *
* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
* for more details. (http://www.gnu.org/licenses/)                           *
* -------------------------------------------------------------------------- */

#include "LinAlgDecl.h"
#include "Matrix.h"
#include "Vector.h"

#include <algorithm>
#include <cmath>

using namespace std;

namespace RSIM
{

/** 
        
Note: This file was obtained from Template Numerical Toolkit (TNT).
Details of TNT can be found on the website http://math.nist.gov/tnt/
This file is minimally modified to make it compatible with RSIM.
The originial file name (available for download on TNT website) 
is jama_eig.h.

<P>
Computes eigenvalues and eigenvectors of a real (non-complex)
matrix. 
<P>
If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
diagonal and the eigenvector matrix V is orthogonal. That is,
        the diagonal values of D are the eigenvalues, and
V*V' = I, where I is the identity matrix.  The columns of V 
represent the eigenvectors in the sense that A*V = V*D.

<P>
If A is not symmetric, then the eigenvalue matrix D is block diagonal
with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
eigenvalues look like
<pre>

        u + iv     .        .          .      .    .
        .      u - iv     .          .      .    .
        .        .      a + ib       .      .    .
        .        .        .        a - ib   .    .
        .        .        .          .      x    .
        .        .        .          .      .    y
</pre>
        then D looks like
<pre>

        u        v        .          .      .    .
        -v        u        .          .      .    . 
        .        .        a          b      .    .
        .        .       -b          a      .    .
        .        .        .          .      x    .
        .        .        .          .      .    y
</pre>
This keeps V a real matrix in both symmetric and non-symmetric
cases, and A*V = V*D.



<p>
The matrix V may be badly
conditioned, or even singular, so the validity of the equation
A = V*D*inverse(V) depends upon the condition number of V.

<p>
        (Adapted from JAMA, a Java Matrix Library, developed by jointly 
        by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
**/

class Eigenvalue{
        private:
                
                /** Row and column dimension (square matrix).  */
                int n;

                int issymmetric; /* boolean*/

                /** Arrays for internal storage of eigenvalues. */

                Vector d;         /* real part */
                Vector e;         /* img part */

                /** Array for internal storage of eigenvectors. */
                Matrix V;

                /** Array for internal storage of nonsymmetric Hessenberg form.
                @serial internal storage of nonsymmetric Hessenberg form.
                */
                Matrix H;
                
                /** Working storage for nonsymmetric algorithm.
                @serial working storage for nonsymmetric algorithm.
                */
                Vector ort;

                // Symmetric Householder reduction to tridiagonal form.
                void tred2();

                // Symmetric tridiagonal QL algorithm.
                void tql2 ();

                // Nonsymmetric reduction to Hessenberg form.
                void orthes ();

                // Complex scalar division.
                double cdivr, cdivi;
                void cdiv(double xr,double xi,double yr,double yi);

                // Nonsymmetric reduction from Hessenberg to real Schur form.
                void hqr2();
        public:

                /** Check for symmetry, then construct the eigenvalue decomposition
                @param A    Square real (non-complex) matrix */
                Eigenvalue(const Matrix& A);

                /** Return the eigenvector matrix
                @return    V */
                void getV(Matrix& V_){
                        V_ = V;
                        return;
                }

                /** Return the real parts of the eigenvalues
                @return     real(diag(D)) */
                void getRealEigenvalues(Vector& d_){
                        d_ = d;
                        return ;
                }

                /** Return the imaginary parts of the eigenvalues
                in parameter e_.
                @param e_: new matrix with imaginary parts of the eigenvalues. */
                void getImagEigenvalues (Vector& e_) {
                        e_ = e;
                        return;
                }


/** 
        Computes the block diagonal eigenvalue matrix.
If the original matrix A is not symmetric, then the eigenvalue 
        matrix D is block diagonal with the real eigenvalues in 1-by-1 
        blocks and any complex eigenvalues,
a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
eigenvalues look like
<pre>

        u + iv     .        .          .      .    .
        .      u - iv     .          .      .    .
        .        .      a + ib       .      .    .
        .        .        .        a - ib   .    .
        .        .        .          .      x    .
        .        .        .          .      .    y
</pre>
        then D looks like
<pre>

        u        v        .          .      .    .
        -v        u        .          .      .    . 
        .        .        a          b      .    .
        .        .       -b          a      .    .
        .        .        .          .      x    .
        .        .        .          .      .    y
</pre>
This keeps V a real matrix in both symmetric and non-symmetric
cases, and A*V = V*D.

        @param D: upon return, the matrix is filled with the block diagonal 
        eigenvalue matrix.
        
*/
                void getD (Matrix& D);
}; // class Eigenvalue
} //namespace RSIM

#endif
// RSIM_EIG_H
